#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdbool.h>
#include <quadmath.h>

#include "../src/myexp.h"
#include "../src/myexp_impl.h"

#if defined(TEST_EXP_V0)
    #define EXP_FUN     exp
    #define EXP_FNAME  "   exp"
#elif defined(TEST_EXP_V1)
    #define EXP_FUN     exp_v1
    #define EXP_FNAME  "exp_v1"
#elif defined(TEST_EXP_V2)
    #define EXP_FUN     exp_v2
    #define EXP_FNAME   "exp_v2"
#elif defined(TEST_EXP_V3)
    #define EXP_FUN     exp_v3
    #define EXP_FNAME  "exp_v3"
#else
   #error must define TEST_EXP_V0 or TEST_EXP_V1 or TEST_EXP_V2 or TEST_EXP_V3
#endif

#define ARY_SIZE(a) sizeof(a)/sizeof(a[0])

void test_exp_func(double *p,int len)
{
  double x, x1,y0,y1;

  for (int i = 0; i < len; i++ )
  {
    x = p[i];
    x1 = x * INV_LN2;
    y0 =(double)expq((__float128)x);
    y1 = EXP_FUN(x);

    printf("exp(%.16lf)=2^%.16lf=%.13la:ref_exp\n", x, x1, y0);
    printf("exp(%.16lf)=2^%.16lf=%.13la:%s\n",x, x1, y1,EXP_FNAME);
  }
}

void test_exp_case1()
{
  printf("\nrun %s()\n", __func__);
  char buf[64];
  double xs[] =
  {
    -1.0/0.0 , 1.0/0.0 /* nan*/,
    -__builtin_inff(), __builtin_inff(),
    -1025.0, -1024.0, 1024.0, 1025.0
  };
  test_exp_func(xs, sizeof(xs)/sizeof(xs[0]));
}

void test_exp_case2()
{
  printf("\nrun %s()\n", __func__);
  int e;
  double f,eps,x,x1,y0,y1;
  //1024/log(2) ~= -1477.31972187
  //2024*log(2) ~= 709.7827128933839968432
  for (int i=-1478; i<=1048;i++)
  {
    f=(double)i*LN2;
    frexp(f,&e);
    eps = ldexp(0.5,e-52);
    for (int j=-1;j<=1;j++)
    {
      x  = f + ((double)j * eps);
      x1 = x * INV_LN2;
      y0 =(double)expq((__float128)x);
      y1 = EXP_FUN(x);
      printf("exp(%.16lf)=2^%.16lf=%.13la:ref_exp\n", x, x1, y0);
      printf("exp(%.16lf)=2^%.16lf=%.13la:%s\n",x, x1, y1,EXP_FNAME);
    }
  }
}

void test_exp_case3()
{
  printf("\nrun %s()\n", __func__);
  double f,x,x1,y0,y1;

  for (f=-10.0; f<=10; f= f +0.01 )
  {
    x  = f;
    x1 = x * INV_LN2;
    y0 =(double)expq((__float128)x);
    y1 = EXP_FUN(x);
    printf("exp(%.16lf)=2^%.16lf=%.13la:ref_exp\n", x, x1, y0);
    printf("exp(%.16lf)=2^%.16lf=%.13la:%s\n",x, x1, y1,EXP_FNAME);
  }
}

void debug_exp()
{
  double x = -1022;
  double f =x*LN2;
  printf("f=%.16lf\n",f);

  double y0 = exp(f);
  double y1 = EXP_FUN(f);
  printf("ref_exp(%.16lf)=%.13la\n",f,y0);
  printf("my_exp(%.16lf)=%.13la\n",f,y1);

  double r0=y0/0x1p-1022;
  double r1=y1/0x1p-1022;
  printf("r0=%.16lf\n",r0);
  printf("r1=%.16lf\n\n",r1);

  f=-42.28197801415663635;
  // f=-42.2819780141566355;
  y0 = exp(f);
  y1 = EXP_FUN(f);
  printf("y0=%.13la\n",y0);
  printf("y1=%.13la\n",y1);

}

void test_exp_fun()
{
  test_exp_case1();
  test_exp_case2();
  test_exp_case3();
}

int main(int argc, char *argv[])
{

#ifdef _DEBUG_
  //debug_exp();
  return 0;
#endif

  test_exp_fun();

  return 0;
}
